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ABSTRACT.  The  problem  of  constructing  a  surface  from  the  information  provided  by 
the  Marr-l*oggio  theory  of  human  stereo  vision  is  investigated.  It  is  argued  that  not 
only  does  this  theory  provide  explicit  boundary  conditions  at  certain  points  in  the 
image,  but  that  the  imaging  process  also  provides  implicit  conditions  on  all  other 
points  in  the  image.  This  argument  is  used  to  derive  conditions  on  possible 
algorithms  for  computing  the  surface.  Additional  constraining  principles  are  applied 
to  the  problem;  specifically  that  the  process  be  performable  by  a  local-support 
parallel  network.  Some  mathematical  tools,  differential  geometry,  Coons  surface 
patches  and  iterative  methods  of  convergence,  relevant  to  the  problem  of 
constructing  the  surface  are  outlined.  Specific  methods  for  actually  computing  the 
surface  are  examined. 
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1.  Introduction 

In  a  recent  article,  Marr  8c  Poggio  [1977]  set  out  a  computational  theory  of 
human  stereo  vision.  One  consequence  of  this  theory  is  that  the  stereo  algorithm  can 
at  host  determine  disparity  only  along  certain  contours  in  the  image.  Two  important 
questions  thus  become:  Is  it  possible  to  reconstruct  a  surface  which  is  consistent 
with  the  stereo  information,  and  if  so,  how  does  one  reconstruct  such  a  surface? 
This  paper  addresses  these  questions  and  discusses  several  pieces  of  mathematics 
relevant  to  the  interpolation  of  surfaces  under  such  conditions. 

The  original  motivation  for  the  problem  lies  in  the  theory  of  human  stereo 
vision:  thus,  we  seek  an  algorithm  which  is  plausible  as  a  human  model,  although 
we  will  not  claim  that  the  resulting  algorithm  is  an  exact  model  of  a  module  of  the 
human  system. 

In  designing  a  plausible  algorithm,  a  number  of  constraining  principles  are 
applied  to  the  problem,  and  are  accepted  as  given.  Any  visual  system  must  be  able  to 
process  largo  amounts  of  input  data;  for  example,  in  the  human  system,  the  central 
part  of  the  visual  field  may  be  considered  as  an  image  which  is  a  thousand  pixels  on 
a  side.  '  At  the  same  time,  it  is  desired  that  the  system  perform  its  computations  in 
real  time;  i.o.  that  it  take  only  a  short  period  of  time  after  receiving  the  input  to 
complete  its  computation.  Since  each  processor  must  take  some  non-negligible 
amount  of  time  to  perform  each  action,  this  essentially  implies  the  use  of 
computations  which  can  be  implemented  in  a  parallel  manner,  using  a  large  number 
of  interconnected  processors.  This  is  the  first  constraint. 

A  second  physical  constraint  on  our  process  is  that  the  "hardware"  which 
implements  it  must  fit  into  a  finite  and  constrained  amount  of  space  (such  as  the 
cranium).  The  use  of  parallel  networks  of  interconnected  processors,  combined  with 
this  constraint  of  physical  space  available,  requires  that  each  processor  not  be 
connected  (o  all  others.  Mather,  there  should  only  be  local  connections  between  the 
processors.  Here,  local  means  not  only  that  the  number  of  connections  be  small,  but 
that  since  we  are  processing  information  whose  underlying  coordinate  system  is  a 
two-dimensional  plane,  the  connections  should  also  be  local  in  a  spatial  sense.  If  the 
support  of  a  function,  defined  on  a  two-dimensional  grid,  is  the  set  of  points  on  the 
grid  which  contribute  in  a  non-trivial  manner  to  the  computation  of  the  function. 


1.  IIm  M/.<  o|  dir  nnliMilii.il  rncpiois  in  the  center  of  (ho  focral  region  is  r  our  lily  20  to  ,15  second* 
ol  -ii'  l  I  V  (.nmsMirt,  Visual  lYrrrpimn  { 1 070].  p.  156).  tlvrr  llic  central  6  drprrrs  of  the  visual  field, 
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then  our  requirement  is  that  the  processors  implementing  our  computation  must 
have  local  (topological  or  metric)  support. 

A  third  requirement,  though  not  as  strong  as  the  first  two,  concerns  the 
complexity  of  the  individual  processors.  We  do  not  want  to  achieve  locality  at  the 
expense  of  computation  time.  That  is,  if  the  computation  can  be  made  local  only  at 
the  expense  of  requiring  each  processor  to  compute  some  complex  function  requiring 
a  long  time  interval  to  complete,  then  there  is  something  fundamentally  wrong  with 
the  computation  as  a  part  of  the  image  processing  system.  Note  that  the  global 
computation  performed  by  the  parallel  network  need  not  be  simple,  it  is  only 
necessary  that  the  individual  processors  not  be  complex.  Connected  with  this  desire 
for  simplicity  is  a  second  desire  for  uniformity,  that  is,  if  possible,  the  individual 
processors  in  the  network  should  be  identical.  However,  this  is  not  as  critical  a 
requirement  as  that  of  parallelism  and  local  support. 

Although  the  original  motivation  for  such  constraints  on  the  algorithm 
arise  from  consideration  of  the  human  visual  system,  they  could  apply  equally  well 
to  other  typos  of  image  processing  systems,  and  are  taken  as  general  constraints  on 
the  computation  wo  are  about  to  investigate.  As  a  consequence,  this  paper  will  not 
suggest  a  particular  algorithm  as  a  model  of  the  human  system.  Rather,  a  number  of 
alternatives  will  be  suggested  and  examined  in  terms  of  their  acceptability  vis-a-vis 
the  algorithmic  constraints  outlined  above. 

The  problem  is  approached  in  the  following  manner.  We  first  determine 
exactly  what  information  is  available  from  the  stereo  algorithm.  It  is  shown  in  the 
second  section  that  the  stereo  algorithm  gives  implicit  as  well  as  explicit  information 
about  the  shape  of  the  surfaces  in  a  scene.  Next,  we  turn  this  information  into 
conditions  and  constraints  on  the  computation.  In  section  three  we  outline  general 
methods  for  satisfying  these  constraints.  This  includes  an  outline  of  some 
differential  geometry  relevant  to  the  description  and  reconstruction  of  surfaces. 
Also,  several  extremal  conditions  for  the  interpolation  of  a  surface  are  suggested. 
.Section  four  outlines  the  Coons  surface  patch  method  for  creating  smooth  fair 
surfaces  from  boundary  conditions.  Finally,  we  design  an  actual  algorithm  which 
satisfies  the  algorithmic  constraints  of  simple,  parallel,  local-support  processes  and 
which  solves  the  problem  subject  to  the  computational  constraints  developed  in  the 
previous  sections.  Section  five  outlines  methods  for  performing  the  computation. 
Section  six  combines  the  previous  sections  and  illustrates  how  to  actually  construct 
the  surface. 
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2.  Information  from  Stereo 

Accord  inf,  to  the  current  computational  theory  of  human  stereo  vision 
( Marr  &  Poggio  [1977]),  the  human  visual  processor  solves  the  stereoscopic  matching 
problem  by  means  of  an  algorithm  that  consists  of  five  main  steps:  (1)  The  left  and 
right  images  are  each  filtered  at  different  orientations  with  bar  masks  of  four  sizes 
that  increase  with  eccentricity;  these  masks  have  a  cross-section  that  is 
approximately  the  difference  of  two  gaussian  functions,  with  space  constants  in  the 
ratio  hl.Vfi.  (2)  Zero-crossings  in  the  filtered  images  are  found,  along  scan-lines 
lying  perpendicular  to  the  orientation  of  the  mask.  Termination  points  of  lines  and 
edges  are  also  localized.  (3)  For  each  mask  size,  matching  takes  place  between  pieces 
of  zero-crossing  contour  of  the  same  sign  and  roughly  the  same  orientation  in  the 
two  images,  for  a  range  of  disparities  up  to  about  the  width  of  the  mask's  central 
region.  Within  this  disparity  range,  Marr  &  Poggio  showed  that  false  targets  pose 
only  a  simple  problem,  (d)  The  output  of  the  wide  masks  can  control  vergence 
movements,  thus  causing  small  masks  to  coinc  into  correspondence.  In  this  way,  the 
matching  process  gradually  moves  from  dealing  with  large  disparities  at  low 
resolution  to  dealing  with  small  disparities  at  high  resolution.  (5)  When  a 
correspondence  is  achieved,  it  is  stored  in  a  dynamic  buffer,  called  the 
?.  1  /2-dimensional  sketch. 

The  .justification  for  this  model  of  stereo  processing  is  well  detailed  in 
Marr  K  Poggio  [  1977],  and  will  not  be  dealt  with  here.  Our  concern  in  this  paper  is 
how  to  use  the  information  provided  by  the  stereo  algorithm  to  construct  a 
representation  of  the  underlying  surfaces  in  the  image.  To  do  this,  one  must 
carefully  consider  what  information  is  actually  provided  by  the  stereo  mechanism. 

The  important  features  of  the  stereo  theory,  from  the  point  of  view  of  this 
pap'-r,  are  the  following.  Each  image  is  convolved  with  a  mask  which  is  a 
two-dimensional  difference  of  gaussians.'  The  convolutions  are  searched  for 
zero-crossings  and  the  position  and  sign  of  such  zero-crossings  are  stored  in  a 
description  of  ihe  image.  These  zero-crossing  contours  form  the  basic  descriptors 
which  are  matched  by  the  stereo  algorithm.  As  a  consequence  the  only  places  in  the 
image  which  may  have  explicit  disparity  information  associated  with  them  are  those 
corresponding  to  zero-crossing  contours. 

I.  In  i  In*  artn.il  iinpli  mi  lit  at  mu  nf  tin  Hum  theory,  a  rirnil.tr  symmetric  difference  of  gaussianft 
in.i-1.  is  iisiiI  i.ithn  than  an  oriented  mask.  The  instil ir.iiion  for  the  use  ol  sueh  masks  is  found  in 
Marr  £  Hildrilli  [1070]. 
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There  is  strong  psychophysical  evidence  for  such  operators  (Marr  &  Poggio 
[  1 977 J,  Wilson  (iiese  [197V],  Wilson  fc  Bergen  [1979]).  What  is  the 
computational  motivation  for  such  operators?  In  the  one-dimensional  case,  such  an 
operator  detects  an  intensity  change;  specifically,  the  points  of  inflection  in 
intensity.  For  the  two-dimensional  case,  we  again  want  to  detect  those  points  which 
correspond  to  a  point  of  inflection,  now  for  some  directional  derivative.  Marr  & 
Hildreth  [1979]  have  shown  that  the  desired  orientation  of  the  directional  derivative 
should  he  such  as  to  coincide  with  the  local  orientation  of  the  underlying  line  of 
zero-cro;  sings.  Under  certain  conditions,  this  orientation  is  the  one  at  which  the 
zero-crossing  lias  a  maximum  slope  and  this  can  be  detected  using  an 
orientation-independent  differential  operator,  the  Laplacian  V?.  Thus,  these 
operators'  essentially  detect  inflections  in  intensity  for  some  resolution  -  the  smaller 
the  mask's  central  excitatory  region,  the  sharper  the  inflection  must  he  in  order  to  be 
detected.  Note  that  these  are  not  discontinuities  in  an  analytic  sense,  but  rather 
these  are  "discontinuities"  with  respect  to  some  resolution  of  the  image. 

From  what  do  such  inflections  in  intensity  arise?  Image  formation  is 
affected  by  four  factors  (Woodham  [  1978]): 

-  imaging  geometry 

-  incident  illumination 

-  surface  photometry 

-  surface  topography. 

Surface  photometry  refers  to  how  light  is  reflected  by  the  object  surface.  It  is 
determined  by  optical  constants  of  the  object  material  and  by  the  surface 
microstructure  (detail  too  fine  to  be  resolved  but  which  causes  observable  effects  in 
the  way  light  is  reflected  at  the  object  surface).  Surface  topography  is  the  surface 
detail  which  is  within  the  resolution  limits  of  the  imaging  hardware.  It  refers  to 
the  gross  object  shape  relative  to  the  viewrer.  If  we  assume  some  illumination  and 
imaging  geometry,  then  inflections  in  the  image  intensities  will  be  caused  either  by 
the  surface  photometry  or  the  surface  topography. 

Thus  sharp  changes  in  color  or  reflectance  of  the  surface,  scratches  in  the 
surface,  sharp  changes  in  the  shape  of  the  surface  all  can  give  rise  to  intensity 
inflections.  These  intensity  inflections  will  cause  zero-crossings  in  the  convolutions 
and  it  is  these  primitive  descriptors  which  are  matched  by  the  stereo  algorithm. 
Given  the  /act  that  inflections  are  caused  by  such  effects  of  the  surface  photometry 


I.  I'lir  del, ills  ol  Mirk  "ril/m  ili'irrior*"  in.iv  In'  fniiinl  hi  Marr  ft  HiMrrlli  [1979]. 
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■'nil  topography,  one  consequence  of  the  stereo  algorithm  is  that  at  best  it  returns 
disparity  values  along  some  set  of  contours  in  the  image.  One  may  explicitly 
determine  depth  or  surface  orientation  only  along  such  contours.  Our  task  is  to 
reconstruct  a  description  of  the  surface  (either  in  depth  or  in  surface  orientation)  at 
all  points  in  the  image. 

In  general,  any  one  of  a  multitude  of  widely  varying  surfaces  could  fit  the 
boundary  conditions  imposed  by  the  stereo  algorithm.  But  to  be  completely 
consistent  with  the  stereo  process,  such  surfaces  must  meet  the  depth  or  surface 
orientation  conditions  along  the  zero-crossing  contours  and  not  give  rise  to  any  other 
zero-crossing  contours  which  do  not  appear  in  the  convolved  image.  This  is  captured 
by  the  following  assertion: 

Assertion:  Places  of  no  information  are  actually  places  of  information. 

By  this,  we  mean  that  for  the  locations  of  the  image  not  associated  with  a 
zero-crossing,  we  may  assume  that  the  underlying  surface  does  not  change  in  a 
radical  way.  An  intuitive  way  of  looking  at  this  claim  is  as  follows.  Suppose  we  are 
given  a  closed  zero-crossing  contour,  within  which  there  are  no  other  zero-crossings. 
An  example  would  be  a  circular  contour,  along  which  the  disparity  is  constant.  One 
surface  which  is  consistent  with  this  set  of  boundary  conditions  is  a  sphere. 
However,  one  could  also  fit  a  highly  convoluted  surface  (for  example  sin(r)/r)  to 
this  sot  of  boundary  values.  Yet  in  principle,  such  a  rapidly  varying  surface  should 
give  rise  to  other  zero-crossings,  since  the  intensity  across  the  surface  should  also 
vary  considerably.  Hence,  we  claim  that  the  set  of  zero-crossing  contours  contains 
implicit  information  about  the  surface  as  well  as  explicit  information,  and  such 
information  can  be  valuable  in  reconstructing  the  surface. 

This  assertion  is  equivalent  to  the  following  statement: 

Assertion:  Except  under  certain  singular  conditions,  a  rapid  change  in  the 
direction  of  curvature  of  a  surface  must  give  rise  to  an  inflection  in 
intensities,  i.c.  you  cannot  hide  a  dip. 

a 

In  order  to  prove  this  assertion,  wo  shall  need  to  first  develop  tools  for 
dealing  with  image  formation.  This  will  allow  us  to  relate  changes  in  the  surface 
orientation  of  a  surface  element  to  changes  in  the  perceived  intensities.  Then  we 
may  state  the  problem  formally  as  a  set  of  lemmas  dealing  with  one  and  two 
dimensional  cases. 


Since  the  intensities  with  which  we  shall  deal  are  caused  by  both  surface 
photometry  and  surface  topography,  it  is  possible  to  have  complex  interactions 
between  the  two  effects.  In  the  case  in  which  both  lactors  have  roughly  equivalent 
effects,  one  could  construct  situations  in  which  the  surface  topography  changes 
radically,  yet  the  surface  photometry  also  changes  sufficiently  so  that  there  are  no 
noticable  changes  in  the  image  intensities.  In  such  situations,  it  seems  unlikely  that 
one  cun  determine  any  information  about  the  shape  of  the  surface  from  the  image 
intensities  themselves.  We  concentrate  instead  on  those  situations  in  which  the 
changes  in  the  surface  topography  dominate  changes  in  the  surface  photometry. 

In  order  to  prove  this  assertion,  we  develop  some  tools  for  dealing  with 
image  formation  [Horn  1970,  1975].  In  the  simplest  case  of  a  single  point  source, 
the  geometry  of  reflection  is  governed  by  three  angles.  The  incident  angle  between 
the  local  normal  of  a  surface  portion  and  the  incident  ray  is  called  i,  the  view  angle 
between  the  local  normal  and  the  emitted  ray  is  called  e,  and  the  phase  angle 
between  the  incident  and  emitted  rays  is  called  g.  The  fraction  of  incident 
illumination  at  a  given  surface  point  reflected  in  the  direction  of  the  viewer  is 
denoted  by  the  reflectance  function  <Mi.e,g).  Most  situations  with  more  complicated 
distributions  of  light  sources  can  be  modelled  by  the  superposition  of  single  point 
sources.' 

We  let  A(x,y,z)  be  the  object  irradiauce  at  the  surface  point  (x.y.z),  scaled 
by  the  ratio  of  image  irradiance  to  scene  radiance.  The  object  irradiance  will  be 
constant  or  obey  some  inverse-square  law  with  respect  to  distance  from  the  source 
for  physical  systems.  The  ratio  of  image  irradiance  to  scene  radiance  is  a  constant 
which  depends  on  the  imaging  system. 

Let  r  =  (x,y,z)  be  a  visible  point  on  an  object,  and  r’  be  the  corresponding 
point  in  the  image.  If  b(r')  is  the  image  irradiance  measured  at  the  image  point  r\ 
then  (Horn  [1970]): 

b(r')  =  A(r)<|>(i,e.g). 

If  we  restrict  our  attention  to  situations  in  which  the  light  source  can  be 
considered  distant  relative  to  the  separation  of  object  and  viewer,  then  the  phase 
angle  is  roughly  constant  and  <£(i,c,g)  can  be  replace  by  the  radiance  function  R(p,q), 
where  p  =  z,  and  q  =  zy  are  the  partial  derivatives  of  the  surface  with  respect  to  the 
two  coordinate  variables,  x  and  y.  It  may  be  that  we  can  in  many  cases  decompose 


1.  for  ,i  more  roMtpIrtr  ilrvrlopinrnl  of  ilir  mallirntat  ms  of  imapr  formation,  srr  Horn  [1970,  1975], 
WooiIImiii  [1978]. 
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the  surf. ice  photometry  from  the  surface  topography  in  such  a  manner,  so  that  the 
image  intensities  depend  only  on  p  and  q.  In  fact,  Woodham  [1978]  observes: 

"  No  matter  how  complex  the  distribution  of  incident  illumination,  for 
most  surfaces,  the  fraction  of  the  incident  light  reflected  in  a  particular 
direction  depends  only  on  the  surface  orientation." 

If  we  also  restrict  our  attention  to  situations  in  which  the  image  projection  is 
orthographic,  then 

b(r')=I(x,y) 

where  I(x,y)  is  the  intensity  value  recorded  in  the  image.  Thus,  the  image  equation 
is  given  by-. 

l(x,y)  =  A(x,y,z)H(p,q). 

Note  that  in  the  case  of  uniform  illumination  and  uniform  photometric  properties  of 
the  surface,  A  is  constant.  However,  in  general,  the  intensities  will  be  a  function  of 
the  position  as  well  as  the  orientation  of  a  surface  patch. 


To  prove  the  assertion,  we  first  examine  the  one  dimensional  case.  We 
wish  to  investigate  the  conditions  under  which  a  bending  of  the  surface  forces  an 
inflection  in  the  intensities.  In  what  follows,  we  assume  that  A,  R  and  z  are  second 
differentiable  functions.  We  begin  with  the  following  lemma: 

Lemma:  Suppose  that  the  surface  and  its  reflective  properties  are 
such  that  changes  in  the  orientation  of  the  surface  dominate  changes  in 
intensity  due  to  changes  in  position  on  the  surface.  If  the  surface  contains  at 
least  two  inflection  points,  then  unless  the  reflectance  R  is  constant  for  the 
values  of  p  involved,  the  intensities  must  contain  an  inflection. 

Proof:  Consider  a  second  differentiable  surface  z(x)  which  contains  at 
least  two  adjacent  inflection  points,  at  x(  and  x2.  Then  z„  and  consequently  p,  are 
zero  at  these  points. 
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Differentiating  the  image  intensity  equation  yields 
I,(x)  =  A„(x)R(p(x))  +  A(x)Rp(p(x))px(x). 

Note  that  if  the  reflectivity  is  completely  independent  of  the  position  in  the  imaf.e, 
then  A  is  constant  and  hence  Ax  -  0,  Then,  since  p,(x()  =  p,(x2)  =  0,  we  see  that 
T,(X|)  =  I,( x?)  =  0.  Since  Ix(x)  has  two  adjacent  zero  points,  the  second 
differentiability  of  the  surface  z  and  the  functions  R  and  A  thus  implies  that  Ix  must 
have  an  extremum  for  some  value  of  x"0,  X|<x0<x2.  In  other  words,  I(x)  must  have 
an  inflection  at  ihis  point.  This  is  provided  I„(x)  is  not  zero  everywhere,  and  this  is 
true  provided  that  Rp  is  not  zero  everywhere,  or  that  R  is  not  constant. 

For  the  more  general  case  of  the  reflectivity  being  a  function  of  the 
position  x,  as  well  as  a  function  of  the  orientation  p,  we  assume  that 
|AX/A|  <<  |Rpp, |  almost  everywhere 

In  other  words,  the  major  cause  of  change  in  intensity  is  changes  in  orientation,  so 
that  the  surface  shape  changes  much  faster  than  the  reflectivity. 

Consider  some  neighbourhood  of  the  point  x(.  Because  p(x()  is  an 
extremum,  there  are  constants  a,.  a2  >  0,  such  that  p(X|-ct|)  =  p(X|+a2)  and  such 
that 

sgn{Ix(x,-rt,))  =  sgn{A(xra|)Rp(p(X|-a,))px(xra|)} 
sgn{I,(x,+«2)}  =  sgn{A(x,+a2)Rp(p(x,+a2))px(x1+<»2)} 

Here,  sgn(x)  is  1  if  x*  is  positive,  -1  if  x  is  negative  and  0  if  x  is  zero. 

Using  the  fact  that 

sgn(xy)  =  sgn(x)sgn(y), 
p(X,-«|)  =  p(x,+«2). 

and  that  A(x)  is  non-negative,  we  see  that 

sgn(Ix(X|-«,))  =  sgn(Ix(X|+«2)) 

if  and  only  if 

sgn(px(x,-«,))  =  sgn(px(x,+a2)). 

But  since  x,  is  an  inflection  point,  this  implies  that 
sgn(p„(X|-a,))  /  sgn(px(x,+«2)). 

Bj-  the  second  differentiability  of  I,  I„  must  be  zero  for  some  point  in  this 
neighbourhood. 

Similarly,  we  may  show  that  I,  is  zero  for  some  point  in  a  neighbourhood 
about  X-..  Arguing  as  before,  we  may  then  show  that  I(x)  contains  an  inflection 
point  except  in  the  case  of  R  being  constant. 
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Lemma:  Suppose  that  the  surface  and  its  reflective  properties  are 
such  that  changes  in  the  orientation  of  the  surface  dominate  changes  in 
intensity  due  to  changes  in  position  on  the  surface.  If  the  surface  contains  one 
inflection  point  at  x,  and  Il(p)  has  an  extremum  for  some  value  p0  such  that 
p(x-)  =  p0.  where  x?^X|,  then  unless  the  reflectance  R  is  constant  for  the 
values  of  p  involved,  the  intensities  must  contain  an  inflection. 

Proof:  The  proof  of  this  lemma  is  very  similar  to  the  previous  one.  As 
before,  we  may  show  that  I,  =  0  for  some  point  in  a  neighbourhood  about  xf. 
Consider  a  neighbourhood  of  the  point  x?.  As  before,  we  may  find  constants  <*|  and 
«2  such  that 

sgn{I,(x2-ft,)}  =  sgn{A(x?-a  ( )Rp(p(x2-a  |  ))p,(x2-a  | )} 
sgn{I„(x2+a2)}  =  sgn{A(x2+a2)Rp(p(x2+a2))p,(x2+a2)}. 

In  this  case,  since  there  is  no  inflection  in  the  surface  here, 
sgn(I,(x2-a, ))  =  sgn(Ix(x2+a2)) 

if  and  only  if 

sgn(Rp(p(x2-oi|)))  =  sgn(Rp(p(x2+o2))). 

Rut  since  R  has  an  extremum  at  p(x2), 

sgn(Rp(p(x2-o,)))  S-  sgn(Rp(p(x2+o2))). 

Arguing  as  before,  wo  see  that  l(x)  must  therefore  contain  an  inflection  point. 

The  only  other  possible  case  of  the  surface  containing  an  inflection  point  is 
if  R  is  nionotonic  and  the  surface  has  only  one  inflection  point.  In  thius  case,  it  is 
not  possible  to  guarantee  that  there  be  an  inflection  in  the  intensities. 

Stated  in  less  formal  terms,  the  above  lemmas  show  that  if  the  major 
changes  in  the  intensities  associated  with  a  particular  surface  are  due  to  the  shape  of 
the  surface  and  not  dun  to  changes  in  the  photometry  of  the  surface,  then  certain 
types  of  changes  in  the  surface  shape  must  give  rise  to  zero-crossings.  This  means 
that  even  if  the  surface  is  painted  with  a  varying  shade  of  color,  so  long  as  the  effect 
of  the  shading  does  not  overwhelm  the  actual  bending  of  the  surface,  there  must  be 
an  inflection  in  the  intensities  corresponding  to  the  change  in  shape  of  the  surface. 
Of  course,  it  may  be  possible  to  paint  a  curved  surface  in  such  a  wa}'  as  to  exactly 
counteract  the  effects  of  the  curving  of  the  surface  on  the  intensities,  thereby 
ensuring  that  no  zero-crossing  will  correspond  to  the  change  in  the  surface  shape. 
We  regard  such  situations  as  special  cases  which  will  be  rare  in  the  real  world,  since 
they  require  a  peculiar  coupling  of  the  surface  topography  and  photometry  which  is 
dependent  on  the  viewer  and  illumination  geometry. 
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The  lemmas  can  bp  extended  to  the  two-dimensional  case.  The  proofs  are 
outlined  below. 

Lemma:  Suppose  that  the  surface  and  its  reflective  properties  are 
such  that  changes  in  the  orientation  of  the  surface  dominate  changes  in 
intensity  duo  lo  changes  in  position  on  the  surface.  Furthermore,  suppose  that 
along  some  direction,  the  surface  undergoes  radical  changes  of  shape,  while 
perpendicular  to  this  direction,  the  change  in  surface  shape  is  negligible.  If 
the  surface  contains  two  inflection  points,  then  unless  the  reflectance  R  is 
constant  for  the  values  of  p  involved,  the  intensities  must  contain  an 
inflection. 

Proof:  The  proof  is  similar  to  that  of  the  one-dimensional  case.  Without 
loss  of  generality,  we  assume  that  the  surface  inflections  occur  along  a  line  parallel 
to  the  x  axis.  In  this  case, 

I,(x,y)  =  A„(x,y)H(p.q)  +  A(x.y)Hp(p.q)p,(x,y)  +  A(x,y)Rq(p,q)q,(x,y). 
Since  the  direction  of  concern  is  parallel  to  the  x  axis,  y  is  constant  along  such  a 
slice  of  the  surface,  say  y  =  y0.  As  before,  we  wish  to  show  that  lx(x,y0)  is  zero  for 
some  point  in  a  neighbourhood  of  the  inflection  point  x,  and  for  some  point  in  a 
neighbourhood  of  the  inflection  point  x2.  Since  we  have  assumed  that  the  effects  of 
P.  dominate  those  of  q„,  we  can  essentially  use  the  same  argument  to  show  that  I, 
changes  sign  within  the  neighbourhood  and  hence  must  be  zero  at  some  point  within 
it. 

Lemma:  Suppose  that  the  surface  and  its  reflective  properties  are 
such  that  changes  in  the  orientation  of  the  surface  dominate  changes  in 
intensity  due  to  changes  in  position  on  the  surface.  Furthermore,  suppose  that 
along  some  direction,  the  surface  undergoes  radical  changes  of  shape,  while 
perpendicular  to  this  direction,  the  change  in  surface  shape  is  negligible.  If 
the  surface  contains  one  inflection  point  at  x,  and  R(p)  has  an  extremum  for 
some  value  p0  such  that  p(x2)  =  p0,  where  x2  *  Xj.  then  unless  the  reflectance  R 
is  constant  for  the  values  of  p  involved,  the  intensities  must  contain  an 
inflection. 

The  proof  again  parallels  that  of  the  one-dimensional  case,  with 
modifications  similar  to  those  used  in  the  proof  of  the  above  lemma. 

We  have  shown  that  except  under  certain  singular  conditions,  a  pair  of 
ini  lections  in  the  surface  must  cause  an  inflection  in  intensities.  So  at  some  level  of 
resolution,  the  convolutions  of  the  image  will  detect  such  Inflections  (if  they  are 
large  enough)  as  zero-crossings.  Hence  the  stereo  output  tells  us  not  only  the  shape 
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of  tilt*  surface  along  (ho  zero-crossings,  but  also  that  the  surface  cannot  curve 
rapidly  or  drastically  between  zero-crossings.  Otherwise,  an  inflection  would  cause 
<i  zero-crossing  and  none  are  evident. 

Thus,  we  claim  that  in  general,  places  which  do  not  have  a  zero-crossing 
contour  must  be  "well-behaved"  in  the  sense  of  curving  as  little  as  possible  and  of 
having  the  sign  of  the  curvature  change  as  little  as  possible.  In  the  next  section,  we 
make  precise  the  notion  of  curvature.  This  will  allow  us  to  use  our  assumption 
about  places  without  zero-crossings  to  design  algorithms  for  reconstructing  a 
surface. 
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3.  Differential  Geometry 

In  order  to  develop  methods  for  reconstructing  the  surface,  we  need  to 
relate  the  shape  of  the  surface  to  inherent  properties  of  the  surface.  More 
importantly,  we  need  to  relate  the  shape  of  the  surface  to  conditions  on  the  surface 
orientation  at  each  point  of  the  surface  or  to  conditions  on  the  depth  values  at  each 
point  of  the  surface.  To  do  this,  we  first  review  some  relevant  pieces  of  differential 
geometry.  For  the  most  part,  these  are  taken  from  Weatherburn  [1927], 

We  have  already  indicated  that  the  surface  (or  piece  thereof)  with  which 
we  are  dealing  has  no  discontinuities  in  depth  and  is  smooth.  This  suggests  that  the 
reconstruction  of  the  piece  of  surface  is  intimately  related  to  the  curvature  of  the 
surface.  We  now  make  the  notion  of  the  curvature  of  a  surface  more  precise. 

3.1  Curvature  of  Curves 

For  a  curve,  the  curvature  at  a  point,  k,  is  defined  in  the  following 

manner, 

A  curve  is  the  locus  of  a  point  whose  position  vector  r  (relative  to  a  fixed 
origin)  is  a  function  of  one  parameter.  In  particular,  that  parameter  can  be  taken  to 
be  the  arc  length  of  the  curve,  s,  measured  from  a  fixed  point  on  the  curve. 

The  (unit)  tangent  to  the  curve  at  a  point  is  defined  by 
t  =  dr/ds  . 

In  the  case  where  the  parameter  of  the  curve  is  taken  to  be  arc  length,  the  tangent 
vector  given  above  is  automatically  a  unit  vector. 

The  curvature  at  any  point  on  the  curve  is  then  defined  as  the  arc-rate  of 
rotation  of  the  tangent.  It  is  also  given  more  formally  by  the  relation 
dt/ds  =  *n 

whore  n  is  a  unit  vector  perpendicular  to  t  and  lies  in  the  plane  spanned  by  the 
tangents  at  that  point  and  a  consecutive  point.  In  other  words,  let  P  be  the  point 
associated  with  the  arc  length  value  s,  and  ]pt  P(  be  the  point  associated  with  the  arc 
length  value  s+ds.  As  ds  toads  to  zero,  the  tangents  at  the  points  P  and  P(  will  define 
a  piano,  such  that  both  tangonts  lie  in  this  plane.  The  vector  n  lies  in  such  a  plane, 
and  is  perpendicular  to  t.  This  vector,  n,  is  the  principal  normal  of  the  curve  at  that 
point,  and  the  plane  containing  two  consecutive  tangents  at  P,  is  called  the  osculating 
plane. 
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An  alternative  method  for  defining  the  curvature  of  a  curve  at  a  point  is  as 
follows.  The  circle  of  curvature  at  a  point  P  on  the  curve  is  the  circle  passing 
through  three  points  on  the  curve  which  in  the  limit  coincide  at  the  point  P.  In 
other  words,  consider  a  point  P  on  the  curve  specified  by  some  value  of  the  arc 
length,  say  s0.  Lot  P,  be  the  point  on  the  curve  associated  with  the  arc  length  value 
s0+«  and  let  p?  he  the  point  associated  with  the  arc  length  value  So-«.  For  any  «  these 
three  points  define  a  circle  passing  through  all  three  points.  As  we  let «  tend  to  zero, 
the  associated  circle  converges  to  the  circle  of  curvature.  Its  radius  p  is  related  to 
the  curvature  of  the  curve  at  that  point  by 
p  =  1  /k  . 

Expressed  in  intuitive  terms,  the  curvature  at  a  point  measures  how 
quickly  the  curve  deviates  from  a  straight  line. 

The  particular  normal  to  the  curve  at  P  that  is  perpendicular  to  the 
osculating  plane  is  the  binormal.  Being  perpendicular  to  the  osculating  plane  means 
that  it  is  perpendicular  to  both  t  and  n,  and  hence  is  parallel  to  txn.  This  unit  vector 
is  denoted  b. 

Just  as  the  curvature  *  measures  the  arc  rate  of  turning  for  the  unit  vector 
t,  the  torsion  r  measures  the  arc  rate  of  turning  for  the  unit  vector  b.  It  can  be 
obtained  from  the  relation 
db/ds  =  -t n  . 

In  intuitive  terms,  the  torsion  measures  how  quickly  the  curve  deviates  from  a 
plane. 


3.2  Curvature  of  Surfaces 

Now  we  can  turn  to  surfaces  and  again  define  the  notion  of  curvature. 

A  surface  is  the  locus  of  a  point  whose  coordinates  are  functions  of  two 
independent  parameters,  u  and  v.  Thus  the  parametric  equations  for  a  surface, 
defined  in  a  Cartesian  coordinate  system  are 

x=f,(u,v)  y=f2(u,v)  z=f3(u,v) 

and  the  surface  is  defined  by  some  function  F(u,v)=0. 

Consider  any  curve  drawn  on  the  surface.  Again,  let  s  be  the  arc  length  of 
the  curve,  measured  from  some  fixed  point  on  the  curve  to  the  current  point  {x,y,z). 
Then  the  tangent  to  the  curve  is  the  vector  {x'.y'.z’),  where  the  ’  refers  to 
differentiation  with  respect  to  s.  Now  the  straight  line  generated  by  the  tangent  to 
any  curve  is  called  a  tangent  line.  In  particular,  all  tangent  lines  at  a  point  {x,y,z) 
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are  perpendicular  to  the  vector  {F„,Fy,Fr}.  (Here  we  use  the  notation  F,  to  denote  the 
partial  derivative  of  F  with  respect  to  x.)  To  see  this,  note  that  F  has  the  same  value 
•it  .ill  points  of  the  surface,  hence  it  remains  constant  along  any  curve  as  s  varies. 
Thus 

F,  dx/ds  +  Fy  dy/ds  +  F,  dz/ds  =  0  . 

Thus  { x'.y’.z'}  and  {F,,l'v,Fr}  are  perpendicular.  So  all  tangent  lines  at  a  point  are 
perpendicular  to  this  vector,  and  thus  lie  in  a  plane  through  {x,y,z}  perpendicular  to 
this  vector.  This  is  the  tangent  plane.  The  normal  to  the  plane  at  the  point  of 
contact  is  the  normal  to  the  surface  at  that  point. 

Any  relation  between  the  parameters,  say  f(u,v)=0  represents  a  curve  on 
the  surface  since  then  r  is  a  function  of  only  one  independent  parameter.  In 
particular,  the  parametric  curves  are  those  for  which 
u  =  constant  or  v  =  constant 
Then  if  we  denote 

i'i  =  c)r/cUi  r?  =  dr/civ 

we  have  that  r,  is  a  vector  tangent  to  the  curve  v  =  constant  at  the  point  r. 
Similarly  for  r2. 

Consider  two  neighbouring  points  on  the  surface,  with  position  vectors  r 
and  r+dr,  corresponding  to  the  parameter  values  (u.v)  and  (u+du,v+dv)  respectively. 
Then 

dr  =  r(du  +  r?dv  . 

Since  the  two  points  are  arbitrarily  closely  spaced  points  on  a  curve  passing  through 
them,  the  length  ds  of  the  element  of  arc  .joining  them  is  equal  to  the  actual  distance 
|d  r  |  between  them.  Thus 

ds2  =  r|2du2  +  2r|T2dudv  +  r22dv2  . 

We  df  f  ine 
F.  =  r,2 
F  =  r,r2 
Ct  =  r22  . 

The  se  quantities  are  called  the  fundamental  magnitudes  of  the  first  order,  and, 
together  with  the  following  quantity,  are  of  use  in  computing  characteristics  of  the 
surface. 

H?  =  EG-F? 

By  definition,  the  normal  to  the  surface  at  any  point  is  perpendicular  to 
every  tangent  line  through  that  point.  Hence  it  is  perpendicular  to  both  rt  and  r?. 
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Thus  the  unit  normal  is  given  by 
n  =  r,xr2/Il  . 

In  a  similar  manner  the  second  derivatives  of  r  are  denoted 
1*1 1  =  d2r/du2  r,2  =  d2r/dudv  r22  =  d2r/dv2. 

The  second  order  magnitudes  of  the  surface  are  then  defined  as 
I,  =  nr,, 

M  =  nr,2 
N  =  nr?2 
T2=  LN-M2  . 

We  can  now  formalize  the  curvature  of  the  surface  itself.  We  have  seen 
how  one  may  define  the  curvature  of  a  curve  at  a  point.  Consider  any  plane  which 
intersects  the  surface  at  a  particular  point  P,  and  which  contains  the  normal  to  the 
surface  at  that  point.  The  result  of  such  a  normal  section  is  a  curve,  and  we  may 
evaluate  the  curvature  of  that  curve  at  the  point  P.  However,  there  are  infinitely 
nianj-  planes  through  P  which  contain  the  normal  to  the  surface  at  P.  Can  we 
identify  any  particular  normal  sections? 

Given  a  point  P  on  the  surface,  any  curve  on  the  surface  passing  through 
the  point  will  have  a  tangent  vector  defined  at  the  point.  The  plane  containing  all 
the  tangent  vectors  for  any  curve  passing  through  the  point  is  called  the  tangent 
plane  for  that  point.  Suppose  we  intersect  the  tangent  plane  with  the  surface,  and 
examine  the  rate  at  which  the  surface  deviates  from  the  plane  along  any  particular 
direction.  We  will  find  that  there  are  two  directions  on  the  surface,  at  right  angles 
to  each  other,  such  that  in  one  direction  the  surface  deviates  the  quickest  from  the 
plane,  and  in  the  other  direction  the  surface  deviates  the  slowest.  Both  of  these 
directions  have  the  property  that  the  normal  at  a  consecutive  point  separated  by  an 
infinitesimal  distance  in  either  direction  meets  the  normal  at  P.  This  means  that  the 
curve  for  a  section  along  one  of  these  directions  has  no  torsion,  and  is  subject  to 
curvature  in  only  one  direction.  These  are  the  principal  directions. 

t 

The  values  for  the  curvature  of  the  curves  obtained  by  taking  normal 
sections  along  these  principal  directions  are  extrema,  in  other  words  as  we  change 
thp  direction  of  the  section,  the  curvature  for  a  normal  section  achieves  a  maximum 
at  one  of  the  principal  directions.  It  achieves  a  minimum  at  the  other  principal 
direction,  hot  the  curvatures  of  these  special  sections  be  and  *b  respectively.  In 
the  case  of  a  plane,  the  curvature  of  each  normal  section  is  identical.  In  this  case, 
any  pair  of  perpendicular  directions  may  he  taken  as  the  principal  directions. 
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To  relate  the  principal  curvatures  to  the  previous  definition  of  curvature 
in  terms  of  radii  of  circles,  note  the  following.  A  curve  on  the  surface  such  that  the 
normals  at  consecutive  points  intersect  is  called  a  line  of  curvature.  The  point  of 
intersection  of  consecutive  normals  along  a  line  of  curvature  at  P  is  the  centre  of 
curvature  at  1‘.  Its  distance  front  P  is  a  principal  radius  of  curvature  and  the 
reciprocal  is  a  principal  curvature. 

Thus  at  each  point  there  are  two  principal  curvatures  k„  and  «cb.  These  are 
the  normal  curvatures  of  the  surface  in  the  directions  of  the  lines  of  curvature. 
Given  the  principal  curvatures,  the  curvature  of  the  surface  can  be  described  in  a 
number  of  ways'.  The  first  curvature  of  a  surface  is  defined  by 
J  =  K„  + 

The  second  curvature  of  a  surface  (also  Known  as  the  specific  curvature  or  the 
Gaussian  curvature)  is  defined  by 
K  =  */h. 

These  are  related  to  the  fundamental  magnitudes  by 
J  =  (EN-3FM+GL)/IIZ 
K  =  (LN-M2)/II?  . 

In  intuitive  terms,  the  first  curvature  is  analogous  to  the  curvature  of  a 
curve,  while  the  second  curvature  is  analogous  to  the  torsion  of  a  curve. 

Any  point  on  the  surface  may  be  defined  by  the  value  of  the  Gaussian 
curvature  at  that  point.  Thus,  an  elliptic  point  is  one  where  K  >  0.  In  other  words, 
normal  sections  through  the  point  are  all  convex  or  all  concave,  the  surface  does  not 
intersect  its  tangent  plane  at  this  point.  An  example  is  any  point  of  an  ellipsoid.  A 
parabolic  point  is  one  where  K  =  0.  An  example  is  any  point  of  a  cylinder.  A 
hyperbolic  point  (or  saddle  point)  is  one  where  K  <  0.  In  other  words,  there  are 
both  convex  and  concave  normal  sections,  the  surface  intersects  its  tangent  plane. 
An  example  is  any  point  of  a  hyperboloid  of  one  sheet. 


3.3  Example 

Let  the  surface  be  represented  by  the  vector  r={x,y,z(x,y)}.  Then  the 
derivatives  are 

r  i  =  {l.O.z,} 

rr  =  {(U.77) 

r|Xr2  =  {-z„-zY,l} 

|r,xr?|  =  (1  +z„?+zy?)l/z. 

Since  r,  and  r2  are  tangents  to  the  parametric  curves,  the  normal  to  the  surface  lies 
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perpendicular  to  both  of  them.  Thus  the  unit  normal  to  the  surface  is 
It  =  (1+Z  ,?+Zy?)'l/Z  {-Zj.-Zy.l}. 

Hence 

n=i+z,2 

F  = 

G  =  l+Zy2 

Hz  =  l+z^+z/ 

Similarly 

r,,  =  {0,0.z„} 
r \2  =  {0,0,z(y} 

1  ?2  {0,0,  Zyy} 

L  =  2,,/H 
M  =  zxy/H 
N  =  zyy/H 

T?  =  H  ?  (z„zyy-z,/) 

Thus,  for  calculation  purposes, 

J  =  d/d\  (zx/(l+ZxZ+Zy2)1/Z)  +  d/dy  (zy/(l+Z„Z+ZyZ)l/Z) 

K  =  (l+Z,Z+ZyZ)'Z  (Z^Zyy-Z,/)  . 

This  example  suggests  that  there  may  he  another  useful  representation  of 
J,  to  which  we  now  turn. 

The  divergence  of  a  vector  is  defined  (as  in  Weatherburn)  as 
divF  =  H  zr,  (GF0-FFv)  +  Hzr2(EFv-FFu)  . 

This  can  Vie  used  to  show  that 
div  n  =  -J  . 

In  the  case  of  a  surface  with  parameters  x,  y 

div  n  =  V  n  =  dii|/dx  +  <m2/dy  n={nltn2,n3}  . 

As  well,  if  the  I.aplacian  is  denoted  by  Vz,  then 
2K  =  n  V?n  +  (Vn)z  . 

For  the  case  of  r  =  (x,y,z(x,y)} 

Vzn  =  {Vzn,.Vzn2,Vzn3} 

Vzn,  =  d?n,/c>xz  +  r)zn,/ayz  . 

We  can  now  make  more  precise  the  earlier  suggestion  that  first  curvature 
of  a  surface  is  analogous  to  the  curvature  of  a  curve,  and  second  curvature  is 
analogous  to  the  torsion.  For  a  surface 
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^2r  =  .In 

n  V?n  +  (V-nf  =  2K 


For  a  curve,  one  has  the  analogous  equations,  using  a  one  parameter 
version  of  the  del  operator 
V  =  td/ds  . 

Thus 


div  n  =  V'ii  =  -k 
V?r  =  <cn 

n  V2n  +  (V-nf  =  -r?  . 


Thus  J  is  analogous  to  k  and  2K  is  analogous  to  -r?. 


3.4  Extremal  Conditions 

One  of  the  initial  constraints  applied  to  the  problem  of  conducting  the 
surface  was  based  on  the  manner  in  which  the  zero-crossing  contours  of  the  stereo 
program  were  formed.  In  essence,  it  claimed  that  for  portions  of  the  image  not 
associated  with  a  zero-crossing,  the  underlying  surface  must  change  in  a  "reasonable" 
manner.  What  does  this  imply  about  the  surface  itself? 

Our  assumption  requires  that  between  zero-crossing  contours,  the  surface 
should  change  as  little  as  possible,  and  in  particular  should  not  have  any  inflection 
points  in  depth  that  are  not  necessary.  This  is  since  such  inflections  should  give  rise 
to  other  zero-crossings  which  are  not  indicated  in  the  stereo  output.  At  the  same 
time,  the  boundary  values  along  the  zero-crossing  contours,  around  some  portion  of 
the  surface,  will  impose  a  certain  amount  of  intrinsic  curvature  to  any  surface 
fitting  them.  One  method  for  finding  a  surface  which  will  fit  the  boundary 
conditions  along  the  zero-crossing  contours,  and  yet  will  not  curve  excessively,  is  to 
require  that  the  averagp  curvature  at  any  point  on  the  surface  be  minimal. 

A  theorem  duo  to  Euler  states  that  if  Kn  is  the  curvature  of  the  curve 
obtained  by  taking  a  normal  section  of  the  surface  with  a  plane  whose  orientation 
relative  to  one  of  the  principal  directions  is  0,  then 
k„  =  k3  cosz0  +  Kb  sinz0  . 

By  taking  all  possible  normal  sections  and  integrating,  one  finds  that  the 
average  curvature  at  a  point  is  given  by  J/2.  Hence,  one  possibile  method  for 
obtaining  a  surface  which  fits  u  set  of  boundary  conditions  and  has  a  particular 
extremal  property  related  to  the  curvature  of  the  surface  is: 

1.  Find  the  surface  fitting  tiic  boundary  conditions  which  minimizes 


Oil  Iricnl  ial  fi'oiix'l  i  v 


-  20  - 


Crimson 


//. l2dx  dy  =  J7(V  n)2dv  dy  =  ff(t fa+Kb)zdx  dy. 

A  surface  which  minimize.';  the  above  integral  for  some  set  of  boundary  conditions  is 
a  surface  which  minimizes  the  average  curvature  of  the  surface  at  every  point. 

How  well  does  (his  condition  satisfy  our  constraints?  As  we  shall  see,  it  is 
possible  to  construct  a  computational  scheme  which  computes  a  surface  satisfying 
this  condition  and  yet  is  consistent  with  the  notion  of  a  local  parallel  algorithm.  As 
to  whether  the  surfaces  so  constructed  meet  our  condition  of  "well-behaved" 
surfaces  we  note  the  following.  In  the  case  of  boundary  conditions  which  are 
consistent  with  a  simple  surface,  such  as  a  plane,  a  cylinder,  a  sphere  or  even  an 
ellipsoid,  the  above  condition  will  load  to  the  correct  surface  in  each  case.  This 
certainly  meets  our  requirement  that  the  surface  behave  smoothly  and  not  change 
any  more  than  required.  However,  if  the  surface  determined  by  the  boundary 
conditions  and  the  above  requirement  has  a  hyperbolic  point,  it  is  possible  to 
minimize  J2  at  such  a  point  without  minimizing  the  principle  curvatures.  This  Is 
undesirable  since  it  could  lead  to  surfaces  which  violate  our  requirement  of  no 
extremes  except  at  zero-crossings. 

One  manner  of  altering  the  above  condition  to  handle  hyperbolic  points  as 
well  is  the  following: 

Kind  the  surface  f  itting  the  boundary  conditions  which  minimizes 
//«;,2+« b2*x  d>'  =  -//(»  V2n)dx  dy. 

Here,  it  is  no  longer  possible  to  minimize  the  integrand  without  minimizing  the 
principle  curvatures  as  well.  Thus  we  see m  to  meet  our  analytic  constraint  better 
than  the  previous  suggestion.  However,  finding  a  method  of  computing  the  surface 
is  not  as  easy.  Note  that 

\i?+Kb2  =  (#a+«b)2-?.«9»(b  =  J?-2K  =  -n-V?n  . 

This  gives  an  alternate  form  of  the  integrand  to  be  minimized.  However,  the 
conditions  which  will  minimize  such  an  integral  give  rise  to  a  non-linear  system  of 
equations.  Whereas  it  is  possible  to  devise  a  computational  scheme  which  is 
consistent  with  the  notion  of  a  local  parallel  algorithm,  it  is  difficult  to  prove  that 
such  an  algorithm  will  converge  to  the  correct  surface.  Such  an  algorithm  may  be 
undesirable  for  this  reason. 

We  have  required  that  the  surface  which  fits  the  boundary  conditions 
curve  as  little  as  possible.  We  achieved  this  in  the  first  case  by  requiring  that  the 
average  curvature  at  every  point  be  as  small  as  possible.  An  alternate  method  of 
keeping  the  curvature  of  the  surface  minimal  is  as  follows. 
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Consider  d  small  segment  of  the  surface,  surrounding  a  point  P.  For  each 
point  within  the  segment,  translate  the  unit  normal  at  (hat  point  to  the  origin  of  the 
coordinate  system.  This  results  in  the  inscription  of  some  region(s)  of  the  unit 
sphere,  if  the  surface  is  smooth  in  the  region  of  P,  then  the  inscribed  portion  of  the 
unit  sphere  is  a  single  connected  region.  If  one  considers  the  ratio  of  the  area  of  the 
inscribed  region  to  the  area  of  the  original  region,  and  takes  the  limit  as  the  region 
surrounding  P  tends  to  the  point  itself,  one  obtains  the  Gaussian  curvature  at  P. 
Thus  the  Gaussian  curvature  also  measures  the  amount  of  "bending"  of  the  surface  at 
a  point.  Thus  a  possible  constraint  is: 

3.  Find  the  surface  fitting  the  boundary  condition  which  minimizes 
fj K2dx  dy  =  ./7k/ s/dx  dy  . 

Minimizing  the  above  integral  results  in  a  surface  which  reduces  the  "bending",  as 
measured  by  the  Gaussian  curvature,  to  a  minimum. 

There  may  be  other  extremal  conditions  that  can  be  applied.  The  critical 
point  is  that  they  ensure  that  the  surface  which  satisfies  them  is  consistent  with  the 
constraints  previously  derived.  That  is,  they  must  not  have  additional  inflections 
other  than  at  zero-crossings,  they  must  spread  the  curvature  of  the  surface  in  a 
smooth  manner  over  the  surface,  and  they  must  fit  the  stereo  boundary  conditions 
along  the  zero-crossings  contours. 

Most  of  the  extremal  conditions  listed  above  are  phrased  in  terms  of  the 
surface  normal  at  points  on  the  surface.  Hence  they  are  best  suited  for  constructing 
a  surface  when  the  boundary  conditions  are  specified  in  terms  of  surface  orientation. 
However,  one  can  also  contruct  a  surface  when  the  boundary  conditions  are  specified 
in  terms  of  distance.  The  next  section  examines  this  case. 
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A.  Surface  Patches 

We  have  indicated  that  our  basic  intention  is  to  fit  a  smooth  surface  to  a 
set  of  boundary  conditions.  Since  the  number  of  possible  surfaces  is  infinite  and 
varied,  we  have  restricted  our  attention  to  those  surfaces  which  satisfy  certain 
extremal  conditions  with  respect  to  surface  curvature.  This  section  outlines  a 
particular  method  for  fitting  smooth  surfaces  to  boundary  conditions.  This  method 
is  especially  suited  to  the  task  of  fitting  a  surface  to  boundary  conditions  involving 
depth  information. 

Coons  [1967]  has  developed  a  method  for  piecing  together  patches  of  a 
smooth  surface  in  such  a  way  as  to  ensure  continuity  of  both  the  function  describing 
the  surface  and  its  derivatives  up  to  some  order  r.  The  method  is  described  below. 

Suppose  that  the  surface  is  parameterized  in  terms  of  u,v  so  that  in  a 
Cartesian  coordinate  system  the  coordinates  are  given  by 
x  =  f(u,v) 
y  =  g(u.v) 
z  -  h(u.v)  . 

Utilizing  Coons’  notation,  let  us  denote  any  point  on  the  surface  for  the 
three  dimensional  case  by 

(uv)  =  [f(u.v),  g(u,v).  h(u,v)] 

Furthermore  we  can  scale  u  and  v  so  that  they  range  from  0  to  1.  So  a  surface  patch 
is  a  segment  bounded  by  four  space  curves,  namely  (Ov),  (lv),  (uO),  (ul). 


In  order  to  define  a  smooth  surface  segment  (or  patch)  between  these  boundary 
curves,  we  will  Mend  the  boundary  curves  together  in  a  smooth  manner.  Let  F0  and 
F|  be  two  blending  functions.  Then  the  surface  equation  is  given  by 
(uv)  =  ,c,(iv)F,(u)  +  E,(u,!)Fj(v)  -  E,SIj(i,j)F,(u)FJ(v)  . 

The  F’s  are  restricted  such  that 
F0(0)  =  1  1 '0(  1 )  =  0 

F,(0)  =  0  Fj(  1 )  =  1  . 


I)i(  f  ri rut  i.il  prmnrl  rv 


-  23  - 


Crimson 


This  ensures  that  the  surface  defined  by  the  surface  equation  above 
contains  its  boundaries  and  corner  points.  In  general,  the  functions  F0  and  F(  are 
taken  to  bo  continuous  and  nionotonic,  but  this  is  not  critical.  The  boundary  curves 
should  be  closed  and  continuous. 

By  imposing  the  additional  restrictions 
F0‘(0)  =  0  F0'(  1 )  =  0 

F,'(0)  =  0  F,*(l)  =  0 

the  slope  of  the  surface  across  a  boundary  has  the  form 
(u0)v  =  (00)vF0(u)  +  ( 10)vF,(u)  . 

Thus  the  slope  across  the  boundaries  depends  on  the  two  end  tangent  vectors  across 
the  boundary  and  on  the  blending  functions.  It  is  independent  of  the  actual  shape  of 
the  boundary  curve.  Thus,  any  two  patches  which  share  a  common  boundary  will 
he  continuous  in  slope  across  this  common  boundary  under  the  above  restrictions  on 
the  blending  functions.  Similarly,  if  the  second  order  derivatives  satisfy 
10”(0)  =  0  F0"(  1 )  =  0 

F|"(0)  =  0  F|"(  1)  =  0 

then  the  patches  will  match  in  the  second  derivative  across  common  boundaries. 

Differentiation  indicates  that 

(00)uv  =  (01  )ov  =  (10)^  =  (00)uv  =  0  . 

In  other  words,  the  cross-derivative  or  twist  vectors  at  the  patch  corners  are  all  zero. 

Not  all  the  surfaces  with  which  we  must  deal  will  have  cross-boundary 
slopes  of  the  above  form,  nor  will  they  all  have  zero  twist  vectors  at  the  corners  of 
the  patches.  In  such  cases,  a  correction  surface  should  be  added  to  the  original 
surface  to  account  for  this.  The  function  of  this  correction  surface  is  to  correct  the 
slopes  across  the  patch  boundary  and  the  corner  twist  vectors  without  changing  the 
shape  of  the  boundary  curves.  In  other  words,  the  correction  surface  changes  only 
slope  and  higher  order  conditions.  This  correction  surface  is  defined  by 
(uv)  =  Z,(iv)ur.,(u)  +  E,(uj)vG,(v)  -  Z1E,(ij)uvGl(u)G,(v)  . 

Adding  the  two  surfaces  together,  where  we  denote  the  initial  surface  by 
f(u.v)  and  the  correction  surface  by  g(u,v),  we  see  that  the  slope  across  the  patch  is 
given  by 

(u0)v  =  f(u,0)v  +  g(u.0)v  . 

Hence,  we  choose  g(u.0)v  and  the  other  cross-boundary  slope  functions  in  g(u,v)  such 
that 

g(u,0)v  =  (u0)v  -  (O0)vF0(u)  -  ( 1 0)VF | ( u )  . 


- 
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The  functions  Gp  and  G,  are  the  slope  blending  functions  and  have  the 
following  end  conditions  in  order  to  ensure  that  the  correction  surface  does  not  alter 
the  patch  boundaries  but  has  the  required  cross-boundary  slopes. 

G0(0)  =  0  G0(  1 )  =  0  Gq’CO)  =  1  G0'(  1 )  =  0 

0,(0)  =  0  G,(  1 )  =  0  G|'(0)  =  0  G,*(l)  =  1  . 

In  a  similar  manner,  one  can  apply  higher  derivative  corrections. 


The  entire  surface  equation  can  be  formulated  in  terms  of  tensors 
( II V  )  —  -[-1  Iq(ii)  F,(ll)  Og(ll)  (»  I  (if )] 

'  tOv)  (00)  (01)  (00)v  (01  )v 

no  no)  (id  no)v  (idv 

(0v)u  (00)o  (01)u(00)uv  (01) 

(10)o  (11  )„  (10)ov  (11) 

This  is  the  equation  for  a  slope-matching,  slope-continuous  surface  patch  with 
arbitrary  boundaries  and  arbitrary  slope  across  the  boundaries. 


■  -1  • 

K0(v) 

F,(v) 

C0(v) 

S’  i*4')- 

Mote  that  the  correction  surface  appears  to  be  essential  in  order  to  achieve 
a  "fair"  interpolated  surface.  This  is  because  the  original  surface  equation  had  zero 
twist  vectors,  yet  most  doubly  curved  surfaces  do  not.  Forrest  observes  [  1968]: 

"It  is  said  that  where  a  series  of  Coons  first  canonical  form  patches 
fi.e.  only  the  basic  surface  equation]  are  fitted  to  an  array  of  points  on  a  car 
body  panel,  and  the  panel  is  cut  by  numerical  control,  the  patches  can  be 
distinguished  on  the  panel  by  a  series  of  flattenings  or  local  distortions;  the 
overall  surface  is  smooth  but  clearly  not  fair." 


We  have  not  yet  specified  the  form  of  the  blending  functions,  and  to  this 
'  we  now  turn.  Let 

111,  UjUaUi,] 

be  a  vector  whose  elements  are  a  set  of  linearly  independent  functions  of  u.  Then 
the  Mending  functions  can  be  specified  by 

(F0(u)  F,(u)  G„(u)  0,(11)]  =  [u,  u2  u3  u4J  M  . 

Thus  we  want  to  specify  M  such  that 


UiLo 

u?Lo 

U3Lo 

1 - 

O 

o 

P 

«i  Li 

«2Ll 

lI3lu-l 

Ufllu-l 

du,/du|o,0 

du2/du|().0 

dU3/du|o.0 

dufl/du|u.0 

_du,/duU, 

du2/rlu(0., 

dug/du^.. 

dUfl/du^, 
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Then  we  have 

F0(u)  =  2u3  -  3u2  +  1 
Fj(u)  =  -2u3  +  3u2 
G0(u)  =  u3  -  2u2  +  xi 
G  |  ( xi )  =  xx3  -  xi2  . 


If  both  basis  vectors,  [xi,  u2  n3  xi4]  and  [v,  v2  v3  v4],  are  taken  to  be  cubic 
basis  vectors,  then  the  resultant  patch  will  be  a  bicubic  surface,  which  is  the 
two-dimensional  analog  of  the  cubic  spline.  Of  course,  many  other  blending 
fxinctions  are  possible,  each  one  resulting  in  different  types  of  smooth  surface 
patches. 
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5.  Convergence  Methods 

The  previous  two  sections  have  dealt  with  the  problem  of  how  to 
reconstruct  a  surface  having  particular  properties  from  a  set  of  boundary  conditions. 
The  reconstruction  itself  has  not  been  formulated,  and  we  now  examine  methods  for 
computing  the  actual  value  of  the  surface  at  various  points,  subject  to  the 
reconstruction  methods,  of  course.  Our  point  of  view  in  handling  this  task  will  be 
to  compute  the  value  of  the  surface  with  as  local  a  support  as  possible. 

Most  of  the  extremal  constraints  on  the  curvature  of  a  surface  can  be 
turned  into  a  set  of  linear  difference  equations,  based  on  a  two-dimensional  grid. 
Similarly,  the  Coons  surface  patch  equation  can  be  turned  into  a  set  of  linear 
difference  equations  based  on  a  two-dimensional  grid.  Thus,  we  examine  iterative 
methods  for  solving  a  set  of  linear  equations.  We  look  at  iterative  methods  rather 
than  elimination  methods  because  we  seek  local,  rather  than  global,  computations. 

Given  a  system  of  linear  equations 
Ax  =  c 

which  is  nonsingular,  we  want  to  find  a  sequence  xk  such  that  xk  converges  to  A*'c. 

An  iteration  is  said  to  be  of  degree  r  if  xk  is  a  function  of  A,  c,  . . . x*.,,.  Usually, 

r=  1  so  that 

xk=Fk(A,  c.  xk.,) 

If  Fk  is  independent  of  k  then  this  is  a  stationary  iteration.  If  Fk  is  a  linear  function 
then  the  iteration  is  linear. 


0.1  Linear  Case 


Suppose  that 

Fk(A.c,xk_|)  =  IlfcX,,.  |  +  vk 

where  Hk  is  a  function  of  A  and  c.  The  solution  should  be  invariant  so 
A**c=H|(A‘lc+VjI 
which  implies 

vk  =  Mkc  where  Mk  =  (I-HJA'1  . 

Thus. 


xk  =  Hkxk.|  +  Mkc  and  Hk  +  MkA  =  I  . 

The  error  associated  with  each  iteration  is 

ek  =  xk-A  'c  where  ek  =  Hkck.|. 

Then  the  system  given  above  is  convergent  for  a  given  initial  error  ejj  If  and  only  if 
K*x  converges  to  0  for  any  x,  where  Kk=HkHk.1...H|.  For  stationary  linear  iterations, 
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Hk=H  and  Kk=Jl\  Thu?  we  need  to  check  that  Hkx  converges  to  0.  It  can  be  shown 
that  ll^x  converges  to  0  for  arbitrary  x  if  and  only  if  each  eigenvalue  X,  of  H  is  such 
that  |A,|<1. 

Bearing  this  in  mind,  we  can  now  develop  .particular  methods  of  iteration. 

5.2  Simultaneous  Displacements  (Jacobi,  1845) 

Given  Ax~c  let  A=E+D+F  where  D  is  a  diagonal  matrix,  E  is  a  lower 
triangular  matrix,  F  is  an  upper  triangular  matrix  and  a,/0  for  all  i.  Suppose  we  are 
given  a  trial  solution  x0  and  we  have  found  xk.|,  (k  =  1,2,...),  then  we  solve 
Kxk.|  +  Dxk  +  Fxk  |  =  c 
or  xk  =  llxk_|  +  Me 

where  II  =  -D  '(E+F)  and  M  =  D'1.  In  order  for  this  iterative  process  to  converge,  it 
is  necessary  and  sufficient  that  the  eigenvalues  X,  of  H  be  bounded  in  modulus  by  1. 
Note  that  if  D  =  6l  then  an  eigenvalue  X,  of  H  corresponds  to  an  eigenvalue  (1-X,)i  of 
A. 

One  set  of  sufficient  conditions  is  given  by  the  following  result,  due  to 
Collate  [1950]. 

If  A  has  diagonal  dominance  and  is  not  reducible,  then  the  method  of 
simultaneous  displacements  converges.  Note  that  A  is  diagonally  dominant  if  and 
only  if 

2,.)|a„|  <  a()  for  all  j 

with  strict  inequality  for  at  least  one  j.  A  is  reducible  if  the  set  {1,...N}  is  the  union  . 
or  two  nonempty  sets  S  and  T  such  that  a,=0  for  all  i«S,  j«T.  In  other  words,  there  is 
a  permutation  matrix  II  such  that 
II  All1  =  TA,  01 

[  A?  A3J 

5.3  Gradient  method 

If  A  is  symmetric  and  definite,  we  can  use  the  gradient  method.  Let 
E(x)  =  x’ax  -  ZcTx. 

Since  E(x)  =  (x-A'lc)1A(x-A‘lc)  -  ctA  'c  it  follows  that  E(x)  attains  a  minimum 
value  -c'a  'c  precisely  when  x=A  'c.  Solving  the  system  Ax=c  is  then  equivalent  to 
finding  the  x  which  minimizes  E(x). 


t 


Given  xp . xk.j.  compute  the  gradient  direction  of  E(x)  at  xk.j 

-VK(x)|,.,  |  =  -2(Axk.,-c)  =  ?,rk.| 

where  rk  is  the  residual.  Since  F.(x)  decreases  in  the  direction  of  rk.|,  we  choose 
Xk  =  Xk-l+nk-|rk-l  . 

We  now  need  to  select  «k  One  possibility  is  to  choose  the  optimal  value, 
which  can  be  shown  to  be 

.Ih: 

rk  i 1  Ark.  | 

A  second  possibility  is  to  lot  ak=a  for  all  k.  Then  the  system  becomes 
=  xk-i  +  «(c-Axk.|) 

and  the  error  term  is  ek=xk-A‘'c  such  that 

ok  =  ok_,  +  a(-Aek.|)  =  (J-aA)ok.|  . 

So  we  have  a  stationary  linear  iteration  with  II=I-aA.  Let  n,  be  the  eigenvalues  of  A, 
A,  those  of  11.  Then  A,=  l-a/i,.  Since  A  is  positive  definite,  we  know  that  n,  >  0  for  all 
i.  The  requirement  of  |A,|  <  1  for  convergence  then  becomes 
0  <  a  <  2/max,  n, . 

For  the  fastest  convergence,  we  must  minimize  max,  |AJ,  so  we  choose  a 
such  that  max,  |l-a/i,|  is  a  minimum. 

If  we  know  that  the  p,  lie  in  the  interval  [a,b],  0<a<b<»,  then  |l-a/i,|  is  a 
maximum  at  the  endpoints.  So  the  best  choice  is  1-oa  =  -(1-ab)  which  implies  1/a  = 
(a+h)/2.  In  this  case, 

1 1 -a/ij  <  (b-a)/(b+a)  <  1  . 

The  gradient  method  is  identical  to  the  method  of  simultaneous 
displacements  whenever  A  is  symmetric  and  definite  with  a  scalar  matrix  as  diagonal 

(n=6i). 

5. 4  Richardson's  method 

We  have  seen  that  in  the  gradient  method 
xk  =  xk.|  +  ok.,rk.,  . 

This  suggests  that  the  relaxation  parameter  o  may  be  a  function  of  k. 

Suppose  A  is  positive  definite  and,  as  before,  the  error  term  is 
ck  =  0-«k-|A)ek., 
ek  -  (I-o0)(l-o,A)...(l-ak.,A)c0  , 


Thus  we  h,ivp 
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«k  =  Pk(A)o0 

Where 

>v.(-v)  =  n0i,<k., 

is  a  polynomial  in  x  such  that  rk(0)  =  1.  Furthermore,  the  zeros  of  Fk  are  1/a,.  So 
choosing  the  constants  a,  is  the  same  as  choosing  a  polynomial  Pk  of  degree  k  such 
that  rk(0)=l.  There  are  a  number  of  ways  of  doing  so;  the  following  is  one  of  the 
more  useful. 

If  the  eigenvalues  of  A  are  p,  then  those  of  Pk(A)  are  Pk(ji,)-  Suppose  that 
M,<[a,bJ  where  0<a<b<»e.  got 
°o  =  2)|<,<n  7,x, 

where  the  x,'s  are  the  eigenvectors  of  A,  and  form  a  basis.  Then 

pk  =  2|<.<n7,Pv.(J‘,)x,  • 

This  suggests  that  one  method  for  ensuring  that  the  error  terms  vanish  quickly  is  to 
make  Pk(p,)  small.  It  can  be  shown  [Forsythe  1960,  p227]  that  one  can  minimize 
ma\-f11,ib  |Pk(x)|  by  choosing 

Fk(x)  =  Tk[(b+a-2x)/(b-a)]/Tk[(b+a)/(b-a)] 
where  Tk(y)  is  the  Ohebyshev  polynomial  and  is  given  by 
Tk(y)  =  cosfk  arccos  y)]  . 

Thus  Pk(x)  is  simply  the  Chebyshev  polynomial  adjusted  to  the  interval 
[d,l>]  and  scaled  to  satisfy  Pk(0)=l.  Let  y0  =  (b+a)/(b-a).  Then 
2Tk(y)  =  [y  +  y2-l),'2]k  +  [y  -  (y2-l)|/2]k  . 

Thus,  as  k  tends  to  infinity, 

|rk(x)|  <  2(y0  -  (y02-l)l/2)k 
The  average  rate  of  convergence  is  then  bounded  by 

-  log(Hekll/lle0ll)l/k  <  -1/k  log  2  +  log  (y0  +  (y02  -  l)1'2  , 

To  compare  this  method  to  the  simultaneous  displacements  method,  let 
P  =  max  n,  /  min  where  p,  are  the  eigenvalues  of  A,  which  is  positive  definite. 
Then  the  rate  of  convergence  for  the  Richardson  method  is 
log  [>'o  +  (yo?-Dl/:?]  =  2(a/b)l/2  <  2P'I/2  . 

The  rate  of  convergence  for  the  method  of  simultaneous  displacements  is 
log  y0  =  2a/b  <2P  1  . 

For  a  process  of  degree  1,  we  cannot  in  general  let  k  tend  to  infinity.  If 
we  use  a  process  such  as 

Xk  =  Xk_,  +ak.,r„., 

we  need  to  know  the  value  of  a,.  Since  the  a,  are  the  reciprocals  of  the  zeros  of 
Pk(x),  we  can  fix  k  at  say  K  =  20,  then  determine  a,  from  the  roots  of  the  Chebyshev 
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polynomial  and  compute  X| . xk.  If  the  error  term  eK  is  not  small  enough,  one  can 

use  the  same  values  of  a,  to  make  another  cycle  of  K  steps.  Such  a  process  is  called 
semi-iterative. 


f).r>  Successive  Displacements  (Gauss-Seidel  [1874]) 

The  method  of  successive  displacements  takes  advantage  of  the 
computation  of  the  components  of  the  vector  as  soon  as  they  become  available.  Thus, 
the  process  is 

F.xk  +  Dxk  +  Fxk.|  =  c 

xk  =  -(D+E)_lFxk  ,  +  (D+Er'c  . 

As  with  the  case  of  simultaneous  displacements,  the  following  results 

hold. 

(CollaU  [1980])  If  A  has  diagonal  dominance  and  is  not  reducible  then  the 
method  of  successive  displacements  converges. 

(Heidi  [1949])  If  A  is  symmetric  and  nonsingular  and  a„>0,  then  the 
method  of  successive  displacements  converges  for  all  initial  states  Xq  if  and  only  if  A 
is  positive  definite. 


5. Cl  Overrolaxntion 

One  can  take  the  method  of  successive  displacements  a  step  further.  For 
processes  such  that  a(i  is  nonzero  for  all  i,  consider 
(E+w_lD)xk  +  [F+(l-w'l)D]xk.|  =  c  . 

The  parameter  w  is  a  relaxation  factor.  If  0<w<l,  then  the  system  is 
underrelaxed.  If  1<m  <2  then  the  system  is  overrelaxed.  In  particular,  if  A  is 
symmetric,  the  diagonal  elements  are  positive  and  the  off-diagonal  elements  are 
non-positive,  then  the  following  is  true. 

If  the  method  of  succesive  displacements  converges  then  the  method  of 
successive  overrelaxation  or  underrelaxation  converges  for  all  <■>  such  that  0<«><2.  If 
the  method  of  successive  relaxation  converges  for  any  0<m<2,  then  the  method  of 
simultaneous  displacements  converges. 

The  value  to  he  assigned  to  the  relaxation  factor  a>  can  be  determined  as 
follows.  Rescale  the  matrix  A  such  that  the  diagonal  elements  are  1.  Let  A  =  -E+I-F. 
Let  p,  be  Ihe  eigenvalues  of  A  and  let  X,  be  the  eigenvalues  of  the  method  of 
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simultaneous  displacements  associated  with  A.  Then  A,  =  1  -  pt.  Let  p,<w)  he  the 
dominant  eigenvalue  of  the  successive  ovorrelaxation  method.  Then  if 

<Ob=  «/[  1  +(J.X|2)|/2] 

tJie  following  relationship  holds  (Kalian  [1958j). 
wb-l  <  h(w)|  <  (wb-l)1/2  . 

Although  the  optimal  relaxation  factor  may  not  be  known,  wb  may  be  close  to  it, 
especially  when  wh-2. 

Although  any  of  the  above  methods  may  be  used  to  solve  a  system  of  linear 
equations,  there  are  various  tradeoffs  associated  with  each  of  them.  For  example, 
although  the  method  of  simultaneous  displacements  tends  to  have  slow  convergence, 
it  can  often  he  implemented  by  a  network  of  processors  which  are  very  local  in 
support.  On  the  other  hand,  the  same  system  of  equations  could  be  solved  by  the 
method  of  successive  displacements.  Here,  the  convergence  is  much  faster,  but  the 
support  of  the  individual  processors  is  much  more  global. 
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t>.  Hoconstruction  of  Surfaces 

Previously,  a  number  of  possible  extremal  conditions  were  suggested  for 
the  construction  of  a  surface.  As  well,  various  iterative  convergence  methods  were 
suggested  for  actually  computing  the  surface.  We  now  consider  the  effects  of 
combining  these  two  factors. 

Consider  the  case  of  finding  the  surface  which  fits  the  boundary 
conditions  and  minimizes 

ffj7  dx  dy  =  fI(K„+Kb)7  dx  dy  . 

Note  that  this  criterion  is  well  suited  for  the  case  of  interpolating  a 
surface  from  surface  normal  information,  since  any  constraints  derived  from  such  a 
condition  directly  relate  to  the  value  of  the  surface  normal  on  the  surface. 

Applying  the  calculus  of  variations  to  this  integral  equation  results  in  the 
necessary  conditions  that  the  partial  derivatives  vanish. 

A/Ax  [Aii|/Ax  +  An?/Ayj  =  A/Ay  [An(/Ax  +  An2/Ay]  =  0  . 

These  Euler  conditions  derived  from  the  calculus  of  variations  can  be  transformed 
into  a  set  of  discrete  conditions  and  thus  into  a  set  of  linear  difference  equations. 
This  allows  one  to  use  the  convergence  methods  discussed  above  to  reconstruct  the 
surface. 

Since  the  partial  derivatives  above  are  bivariate,  the  iterative  matrix  of  the 
set  of  linear  difference  equations  is  block  tridi3gonal.  As  a  consequence,  the  matrix 
as  it  stands  can  be  shown  to  be  divergent.  However,  if  we  transform  the  system  Into 
a  one-dimensional  system  by  means  of  some  simple  matrix  manipulation,  the  system 
becomes  convergent,  and  moreover,  is  very  local  in  nature.  Thus,  the  system  for 
minimizing  the  integral  of  Jz  satisfies  our  condition  of  having  computations  which 
are  local  in  .support.  Moreover,  the  resulting  process  could  easily  be  implemented  in 
a  parallel  network  of  local  support  processors. 

Any  surface  constructed  under  such  a  scheme  can  be  shown  to  be  locally 
composed  either  of  planes,  cylinders,  spheres  or  ellipsoids.  Thus  the  only  basic  type 
of  surface  not  handled  by  this  method  is  hyperboloids. 

It  was  suggested  earlier  that  one  method  for  overcoming  this  handicap 
would  he  to  minimize  the  integral  of  k 7  +  k7  since  such  a  surface  could  contain 
hyperbolic  points.  When  one  applies  the  calculus  of  variations  to  such  a  system,  a 
nonlinear  set  oi  equations  are  obtained,  and  the  convergence  methods  described  above 
are  no  longer  applicable. 
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Thus,  given  the  constraint  of  requiring  a  computational  system  which  is 
local  in  support  and  can  be  computed  by  iterative  methods,  a  good  candidate  for 
reconstructing  surfaces  which  fit  a  certain  set  of  boundary  conditions  and  elsewhere 
contain  no  inflection  points  is  to  find  the  surface  minimizing 
//•  l?dx  ily  =  //(V  n)?dx  dy  =  //(«0+Kl))?dx  dy. 


As  a  second  example,  let  us  consider  how  to  fit  a  smooth  surface  to  a  set  of 
boundary  conditions,  using  the  Coons  surface  patch  method.  For  simplicity,  the 
boundary  is  taken  to  be  a  square  along  which  the  distance  to  the  surface  is  known 
up  to  a  scale  factor.  The  interior  of  the  square  is  of  size  n.  The  basis  vectors  for  the 
Coons  patch  are  cubics  and  the  underlying  parameters  are  taken  to  be  the  axis 
variables,  x  and  y. 

We  treat  the  underlying  system  as  a  two-dimensional  grid  so  that  we  must 
compute  the  depth  of  the  surface  at  each  point  on  an  n*n  grid.  Suppose  we  know  the 
values  of  the  surface  at  the  corners  of  a  particular  patch,  and  we  wish  to  compute 
the  value  for  the  center  of  the  patch,  that  is  for  X  -  y  -  1/2.  Then 
F0(  1/2)  =  Fj(  1/2)  =  1/2  and  Gq(  1/2)  =  -G,(l/2)  =  1/8. 


The  derivatives  can  he  approximated  by  discrete  differences  between 
points  on  the  grid.  By  evaluating  the  tensor  form  of  the  surface  equation,  we  find 
that  the  value  of  the  midpoint  is  given  by 

(.5  .5)  =  .5{(0  .5)  +  (1  .5)  +  (.5  0)  +  (.5  1)} 

-,2f»{(0  0)  +  (0  1)  +  (1  0)  +  (1  1))  . 


We  can  treat  each  surface  patch  as  that  patch  defined  by  a  3*3  piece  of  the 
grid.  Then,  combining  the  surface  equation  for  each  point  of  the  grid,  we  obtain  the 
linear  system 

Ax  =  c 


where  the  constant  vector  c  is  obtained  by  applying  the  above  equation  to  the 
boundary  points  for  the  case  of  a  point  next  to  the  boundary.  Here  the  matrix  A  has 
a  tridiagonal  block  form 


A  = 


A, 

-,6A| 


-.5A| 

A,  -,SA| 


-.6A|  A, 
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Applying  the  Jacobi  method  of  simultaneous  displacements  to  this  block 
system,  yields  the  system 

xk.,  =  Bxk  +  D''c 

with 


In  other  words,  at  each  stage  of  the  iteration,  each  point  in  the  grid 
assumes  the  average  of  its  neighbours  in  the  y  direction,  plus  some  constant  term 
which  reflects  the  influence  of  the  boundaries  in  the  x  direction.  The  constant  term 
is  given  by  D"'c  where  Af1  =  [a„]  with 

a()  =  2/(n+l)  i  (n+l-j)  i<j 

=  2/(n+l)  j  (n+l-i)  i>j 

Provided  the  boundary  of  the  entire  piece  of  the  surface  is  closed,  i.e. 
fixed  depth  values  are  known  along  a  closed  contour  surrounding  the  piece  of  the 
surface,  the  matrix  B  satisfies  the  Chcbyshev  recurrence  relation  and  hence  has 

eigenvalues  oos(kir/n+l)  for  k  =  1.2 . n.  Thus  the  convergence  rate  for  this  process 

is  cos(7r/n+l).  and  the  asymptotic  rate  is  roughly  n'2,  which  is  slow.  However,  notd 
that  the  support  of  each  computation  is  small,  since  each  point  of  the  grid  obtains  its 
value  by  examining  the  values  of  only  two  neighbours. 
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7.  Summary 

This  paper  has  been  a  discussion  of  several  pieces  of  mathematics  relevant 
to  the  interpolation  of  a  surface  from  a  set  of  boundary  conditions.  The  first  section 
outlined  the  problem  as  posed  by  the  Marr  &  Poggio  theory  of  stereopsis.  It  was 
shown  that  the  stereo  process  imposes  both  explicit  constraints  along  the 
zero-crossing  contours  obtained  by  processing  the  image,  and  implicit  constraints 
elsewhere.  These  implicit  constraints  in  particular  included  the  requirement  that 
the  surface  not  change  its  curvature  in  a  radical  manner  for  locations  not  associated 
with  zero-crossings.  The  second  section  dealt  with  some  relevant  details  of 
differential  geoemetry,  in  particular,  those  aspects  dealing  with  the  curvature  of 
surfaces.  The  third  section  outlined  the  method  of  Coons  for  fitting  fair  surfaces  to 
boundary  conditions.  The  fourth  section  sketched  various  iterative  schemes  for 
computing  the  surface,  and  the  fifth  section  tied  all  of  these  notions  together  and 
sketched  the  reconstruction  method. 
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